#include<iostream>
using namespace std;
#include<cmath>

#include"library.h"

double Jacobi_Iter_Solve_Linear(int n, double **a, double *b, double *x, double precision){

	int i,j;
	double y;
	double change = precision +1;
	double *newx = new double [n];
	for(i=0;i<n;i++)x[i] = 0;// initialization

	int count=0;
	while(change > precision){
		change = 0;
		for(i=0;i<n;i++){
			y=0;
			for(j=0;j<n;j++){
				if(j!=i) y += a[i][j] * x[j];
			}
			newx[i] = (b[i] - y)/a[i][i];
			change += fabs( newx[i] - x[i] );
		}
		for(i=0;i<n;i++) x[i] = newx[i];
		count++;
	}
	cout<<"iteration time = "<<count<<endl;
	delete [] newx;
}

int main(){

	int i,j,n=4;
        double A[4][4] = { {5,1,1,1}, {1,6,2,2}, {1,2,8,3}, {1,2,3,10} };
	double **a = new double *[n];
	for(i=0;i<n;i++){
		a[i] = new double [n];
		for(j=0;j<n;j++){
			a[i][j] = A[i][j];
		}
	}
	double b[4] = {8,11,14,16};

	double *x = new double [n];
	
	Jacobi_Iter_Solve_Linear(n, a, b, x, 1E-9);
	cout<<"solution: ";
	for(i=0;i<n;i++)cout<<x[i]<<","; cout<<endl;

	return 0;
}
